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We investigate the dynamics of a phonon-mediated superconductor driven out of equilibrium. 
The electronic hopping amplitude is ramped down in time, resulting in an increased electronic 
density of states. The dynamics of the coupled electron-phonon model is investigated by solving 
Migdal-Eliashberg equations for the double-time Keldysh Green’s functions. The increase of the 
density of states near the Fermi level leads to an enhancement of superconductivity when the 
system thermalizes to the new state at the same temperature. We provide a time- and momentum- 
resolved view on this thermalization process, and show that it involves fast processes associated 
with single-particle scattering and much slower dynamics associated with the superconducting order 
parameter. The importance of electron-phonon coupling for the rapid enhancement and the efficient 
thermalization of superconductivity is demonstrated, and the results are compared to a BCS time- 
dependent mean-field approximation. 

PACS numbers: 74.90.-l-n, 74.40.Gh, 78.47.J- 

I. INTRODUCTION 

Light control of structural and electronic properties 
of solids is a tantalizing prospect of ultrafast materials 
science^^^. In pump-probe experiments, a short pump 
laser pulse drives a solid out of equilibrium. The ensu¬ 
ing dynamics is monitored with a second probe pulse at 
well-defined delay times. Pump excitations at optical fre¬ 
quencies usually create electron-hole excitations, which 
can be used to study transient dynamics in a variety of 
correlated materials^’^, like Mott or charge-density wave 
insulators^^®, or superconductors®^^®. In contrast, lower 
frequency mid-IR or THz lasers can excite the system 
in resonance with structuraP® or other collective modes. 

In particular, intense THz light pulses enable a mode- 
selective vibrational excitation^®, opening up the field of 
“nonlinear phononics”®’^^^^®. 

A lattice deformation can be induced that lasts for hun¬ 
dreds of femtoseconds^®^^®, which has been suggested as 
a basis for light-enhanced superconducting-like nonequi¬ 
librium states^®^^"^. Thus the important question arises 
how fast the electrons in a solid can follow a nonadi- 
abatic change of the lattice structure. In particular, 
the situation is unclear for slow collective modes in a 
symmetry-broken ordered state, such as a superconduc¬ 
tor or a charge-density wave. 

Theoretically, the order parameter dynamics in 
purely electronic models has been investigated in 
BCS mean-field theories for superconductors®®^'^^, and 
in nonequilibrium dynamical mean-field theory for 
antiferromagnets^^’^®. In contrast to such closed systems. 


where the electronic energy is conserved after the ex¬ 
ternal perturbation, in electron-lattice systems energy is 
transferred between electrons and phonons via electron- 
phonon (el-ph) coupling^^. The electronic relaxation in 
electron-phonon models has been theoretically investi¬ 
gated using a variety of methods"^®^®®. 



FIG. 1. THz pump enhances superconductivity, (a) 

Through a lattice distortion, the electronic hopping amplitude 
decreases from Jo to Jf (red shaded area). As a consequence, 
the superconducting order parameter Ao is boosted to a larger 
value A(t) > Aq. At longer time scales (blue shaded area) the 
order parameter approaches its thermal value A/ correspond¬ 
ing to J/ in the presence of efficient electron-phonon (el-ph) 
coupling, (b) Sketch of equilibrium order parameters Aq and 
A/ corresponding to Jo and J/ < Jo, respectively, leading to 
a larger critical temperature Tcj > Tc,o. 

In this work we investigate the nonequilibrium dynam¬ 
ics of a phonon-mediated superconductor induced by a 
transiently modified electronic structure through nonlin¬ 
ear phonon coupling. We consider a tight-binding el-ph 
Hamiltonian which contains both a retarded pairing in- 
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teraction mediated by phonons as well as dissipation of 
heat into the lattice. The light-induced lattice distortion 
is accounted for by a change of the electronic hopping 
amplitude Jq to a smaller value J/ on a typical time 
scale of fractions of a picosecond. Due to this change the 
electronic density of states close to the Fermi surface is 
enhanced, which results in an increased equilibrium order 
parameter Af (see Fig. 1) in the weak-coupling regime 
assumed throughout this work. 

Out of equilibrium, the order parameter is therefore ex¬ 
pected to increase if the change is slow enough and not 
too much energy is deposited into the electronic degrees 
of freedom. Since typically the time scale of the lattice 
distortion - while much longer than the bare electronic 
time scale - is rapid compared to the slow collective dy¬ 
namics of the superconducting condensate, understand¬ 
ing the response of the superconducting order parameter 
A(t) to such a relatively fast change is of great impor¬ 
tance. We show that even for this rapid change of the lat¬ 
tice structure, the superconducting order parameter can 
be drastically enhanced. The dynamics can be separated 
into two different regimes: (i) the short time dynamics 
of the order parameter, which can approximately be de¬ 
scribed by BCS theory, and (ii) the intermediate to long 
time dynamics, where el-ph scattering and the relaxation 
of energy into the dissipative phonon bath dominate. Im¬ 
portantly, the phonon dissipative channel is essential for 
asymptotically reaching the final thermal value. Surpris¬ 
ingly, very fast nonadiabatic ramps are predicted to lead 
to quick enhancement of superconductivity on very short 
time scales in the presence of dissipation. 

The paper is organized as follows: Section II contains 
model and methods. In Section III the main results are 
presented. These results are put into context with con¬ 
clusions and an outlook in Section IV. The Appendix 
contains more detailed information about the BCS for¬ 
malism at finite temperature and additional results on 
the intermetiate time regime. 


II. MODEL AND METHODS 
A. Electron-phonon Hamiltonian 

We investigate the electron-phonon Hamiltonian 
'H = Y^ e{k, t)cl^Ck^ + 

k<7 q,J 

~ (^<3,7 + ^^,7) (^) 

<?,7.o- 

with fermionic creation operators for dimension¬ 
less momentum k = {kx^ky) and spin cr = t)-!- on a 
two-dimensional square lattice with dispersion e{k,t) = 
—2J{t){coskx + cosky) at half-filling. This choice of 
band filling is made for numerical convenience. Away 
from particle-hole symmetric filling, the chemical poten¬ 


tial would have to be adjusted to keep the filling fixed at 
different temperatures, which we avoid. 

The time dependence of the electronic hopping ampli¬ 
tude J{t) mimics a deformation of the lattice induced 
via a nonlinear coupling to an IR active optical phonon 
driven by the THz light pulse^®’®°. Thus, the excited 
phonon is treated classically. We assume for t < r a 
linear ramp J{t) = Jq + [Jf — Jq)^ and for f > r the 
constant J{t) = Jf with Jq = 0.25 eV, J/ = 0.20 eV, 
and ramp time r. The change of the hopping parame¬ 
ter by 20 % is rather large, but not out of reach for an 
experimental realization^®’^^’®^. Furthermore, in an ex¬ 
periment the deformation of the lattice typically lasts for 
several picoseconds, and we focus on the dynamics within 
this time frame. Energies are measured in eV, and time 
scales in fs, using h = 0.658 eVxfs. 

The electrons are coupled to branches ( 7 ) of phonons 
with bosonic creation operators 6 ^.^, energy and 
electron-phonon coupling g^. These quantum phonons 
model the different relaxation channels present in the 
material and should not be confused with the exter¬ 
nally excited phonon mentioned above. We consider 
a dominant optical phonon at Hopt = 0.1 eV, which 
induces superconductivity, and a continuum of acous¬ 
tic low-energy phonons. The distribution of the acous¬ 
tic phonons is given later. We use a reference set of 
electron-phonon couplings, labeled by set (i), and an¬ 
other set with the same spectrum but reduced coupling 
strengths labeled by set (ii). The parameters used for 
the different sets are listed in Table I. We solve this 
model in the Migdal-Eliashberg approximation®^^®'^ with 
a local, self-consistent self-energy for the electrons, and 
treat the phonons as an infinite heat bath at equilib¬ 
rium. The effective phonon spectra weighted by el-ph 
coupling (Eliashberg functions) for case (i) and a param¬ 
eter set without acoustic branch are shown in the inset 
to Fig. 7(a). 

The important difference between BCS mean-field 
treatments, possibly including phenomenological damp¬ 
ing, and the present approach is that we explicitly treat 
a double-time self-energy with nonzero imaginary part, 
which involves true correlation and memory effects and 
accounts for a frequency structure of the phonon spec¬ 
trum. This will be shown to be important in this work in 
the context of absence or presence of low-energy acous¬ 
tic phonons. From the point of view of superconducting 
pairing, the explicit treatment of the retarded self-energy 
also gives a frequency structure to the anomalous self¬ 
energy and thus defines a natural energy cutoff to pairing 
while allowing for a fully gauge-invariant theory. Such a 
cutoff can be introduced in BCS-like theories as well, but 
is often neglected because it violates gauge invariance. 


B. Time evolution 

The time evolution is obtained from solutions of the 
Kadanoff-Baym-Gor’kov equations®^’®® in the Keldysh 
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Green function formalism, described in detail below. We 
choose initial conditions that put the system in the su¬ 
perconducting initial state below T^. We ignore the com¬ 
peting instability towards charge-density wave order at 
half-filling, which is always possible within a mean-field 
scheme. The time-dependent order parameter A(t) is 
defined by 


A(t) ^ Efc/fcW 

Ao "Efc/fc(0) 

using the dimensionless momentum-resolved anomalous 
expectation value fkit) = = {c-ki{t)ck^{t)). 

The initial value Aq = A(t = 0) and final value A/ are 
obtained from the anomalous component of the equilib¬ 
rium self-energy, including energy band renormalization 
with quasiparticle weight Z due to el-ph coupling (see 
below). 

Our choice of a large el-ph coupling A, which results 
in a large value of the order parameter compared to real 
materials®®, is motivated by the times we can reach in 
the numerical simulations. Even though the Migdal- 
Eliashberg approximation is not expected to be quantita¬ 
tively accurate in this regime, the generic effects observed 
should remain valid. 


C. KadanofF-Baym-Gor’kov equations 


We employ the KadanofF-Baym-Gor’kov formalism and 
its application to the superconducting state in the el-ph 
model as described in Ref. 64. We utilize the standard 
two-time Keldysh formalism®®, where the contour Green 
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functions are 2x2 matrices in Nambu space, 






(3) 

(4) 


where t and t' lie on the Keldysh contour, and 7c is the 
contour time-ordering operator. In the following, we use 
units with = /cb = 1. 

The matrix equations of motion with a contour self¬ 
energy ff' to be specified below are 


{idtTo - ek{t))Gk{t,t') = 5^{t,t')To+ 

J^dz (5) 

with 


efc(i) 


€^{k,t) 0 \ 

0 -e 4 ,(-fe,t) J 


( 6 ) 


where fo is the identity matrix. 

On the Keldysh contour, the Langreth rules can be ap¬ 
plied to separate the contour equation into the following 
components: the Matsubara (M), lesser (<), and greater 
(>) Green functions, as well as the mixed real-imaginary 
]/[ Green functions. The various components can be 
transformed or combined into others via the relations 


= (7) 

G^{-iT,ty = G'^{-i{(3 - T),t). (8) 


The equations of motion, letting the contour start at ^min, 
are 


[ - drTo - ek{t^in)]G^^ (t) =iS{r)To - if dz Il^(r - z)G^{z), 

Jo 

pt pj3 

[idtfo - ikit)]G\{t,-iT) = / dz z)G\{z,-iT) - i / dz -iz)G^{z - t), 

'^0 

pt pt' pl3 

[idtfo - ekit)]G^{t,t') = dz z)G^{z,t') + dz z)Gkiz,t') - i dzE^{t, 

J tunin J tm\n JO 


(9a) 

(9b) 


-^z)G[{-iz,t'), 

(9c) 


These equations are solved on the contour by using mas¬ 
sively parallel computation and a time stepping algo¬ 
rithm for integro-differential equations as described in 
Ref. 65 and numerical details are given in Sec. HE. 

Note that we can also drop the spin index f, 7 on the 
normal components of the Nambu Green function, since 


in the absence of magnetic order we have G‘^^{t,t') = 
Gfe 4 ,(i; t') = G^(t, t'). This relation is used for the calcu¬ 
lation of the momentum-resolved normal and anomalous 
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densities 


Ukit) = ( 10 ) 

fk{t) = -iF<{t,t). (11) 


D. Migdal-Eliashberg approximation to 
electron-phonon coupling 


In this work, we employ the Migdal-Eliashberg approx¬ 
imation to electron-phonon coupling of electrons to the 
phononic relaxation channels which are explicitly treated 
in our calculation. These phonons should not be confused 
with the classical phonons involved in the THz driving 
and nonlinear phonon excitation processes. We use a per¬ 
turbative treatment of the electronic self-energy at the 
self-consistent Born level 

= i f dfl a‘^F{n) fa 

( 12 ) 

Here fa is the z Pauli matrix in Nambu space, and 
Gy^^{t,t') = G^{t,t') the local Green function. 

The quantum phonons are kept at hxed equilibrium 
temperature neglecting the phonon self-energy. The 

Keldysh propagator for a single phonon mode at energy 
n is given by 

=-i[nBWn) + 1 - 0c(t,t')]e*^^‘"‘'^ 

-i[nBWn) + (13) 

where UBix) is the Bose function nB{x) = [e^ — 1] 

P = [kBT)~^ the inverse temperature, and is 

the contour Heaviside function. 

The relevant phonon spectrum in Migdal-Eliashberg 
theory is the Eliashberg function^^, defined here for the 
case of a spectrum of local phonon modes 

a^Fin)=Y,\g^\^Sin-n^). (14) 


In practice, we use a model function 

a^F{n) =a^Fopt{n) + (15) 


with an optical branch modelled by a Lorentzian 


a^Foptin) 


2 <5opt 

^°PV((H - Hopt)2 + ^o%) 


(16) 


and an acoustic branch with proper behavior at low 
energy and a cutoff at 2Hacou modelled by 


a'Facou(O) = ^ sin2 
i ^aroii 


ttH 

20 


0(2Oacou - n). 


(17) 


Finally, we estimate a dimensionless electron-phonon 
coupling parameter 


aReEfi(w) 


duj 


ImEff (iwo) 


a ;=0 


Wo 


(18) 


where E)'^ is the normal component of the Matsubara 
self-energy 


rP . _ 

E(iw„) = -i dr e*‘^"^E(T) 

Jo 


with imaginary frequency iujn = i(2n + l)7r//3. 
The renormalized quasiparticle weight 


Z = 


1 

1 + A 


(19) 


( 20 ) 


reflects the effective mass change near Ep induced by el- 
ph coupling. This renormalization is taken into account 
in computing the Bogoliubov dispersions used in Figs. 3 
and 4, and for the initial equilibrium order parameter 

Ao = ZE)^(iwo), ( 21 ) 


where E^'f is the anomalous component of the Matsubara 
self-energy. 

The individual contributions from the phonon modes 
are estimated via 

Aopt = 2N{Ef) dn (22) 

and equivalently for Aacou- Here we use the average elec¬ 
tronic density of states N(Ep) in a window ±Hopt around 
the Fermi level at the initial equilibrium. This expression 
is only strictly valid in the weak coupling limit, where 
the self-energy contributions from optical and acoustic 
branches add up to the total self-energy. 

The detailed parameters used for the runs in this pa¬ 
per are displayed in Table I. We note that the values 
of A used in this work are relatively high in order for 
weak-coupling perturbation theory to be a quantitatively 
accurate approximation. This choice is motivated by nu¬ 
merical feasibility. In order to see a crossover from nona- 
diabatic towards adiabatic behavior, the time scale asso¬ 
ciated with the initial order parameter should not be too 
large compared to the time scales we are able to reach 
in the simulation. A rather large value of A guarantees 
that we have a sizeable order parameter as well as rel¬ 
atively short relaxation times. We expect that using a 
more moderate value of A would change the quantitative 
results, but not the conclusions drawn from our work. 

A few words are in order regarding the choice of the 
heat bath approximation for the phonons, leaving them 
at their initial thermal equilibrium. This can be justified 
by three arguments: (i) The THz-induced lattice modifi¬ 
cation leading to a change of the electronic hopping am¬ 
plitude is a subtle excitation of the electrons, compared 
to the immediate optical excitation of electron-hole pairs 
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Parameter set 

(i) 

(ii) 

w/o acoustic 

Dopt [eV] 

0.100 

0.100 

0.100 

[eV^] 

0.040 

0.032 

0.052 

5opt [eV] 

0.001 

0.001 

0.001 

Aopt 

0.64 

0.51 

0.83 

^acou [sV] 

0.050 

0.050 

- 

flLou [eV^] 

0.025 

0.020 

- 

Aacou 

1.08 

0.86 

- 


TABLE I. Parameters for phonon spectra and electron- 
phonon coupling parameters used in the paper. The param¬ 
eter set without acoustic mode, used for the comparison in 
Fig. 7, has larger coupling strength than set (i) in order to 
match the resulting critical temperature and value of Aq. 


with higher energy photons. Thus, the excess energy in 
the electronic subsystem is relatively small, (ii) We use 
a continuous spectrum of phonons rather than just a sin¬ 
gle sharp mode. The heat capacity of this spectrum of 
phonons is supposedly large, hence the heating of the 
phonons via the transfer of the small amount of excess 
energy is practically negligible, (iii) In reality, if there is 
heating of the phonons, the electrons will thermalize to¬ 
wards the enhanced lattice temperature - at least within 
an effective temperature description - before the excess 
heat is transferred from the irradiated sample surface to 
the bulk crystal. The heating of the lattice strongly de¬ 
pends on the sample and the precise experimental con¬ 
ditions, effects that are beyond the scope of our present 
model study. 


E. Numerical details 

The two-dimensional Brillouin zone was discretized 
with a numerical grid of 80 x 80 momentum points. Cal¬ 
culations were performed on a reduced 1/8 zone with 
a total of 820 momenta. The Kadanoff-Baym-Keldysh 
contour was discretized, with step sizes of St = 0.1 h 
(eV)“^ and St ss 0.21 h (eV)“^ « 0.14 fs for imaginary 
and real times axes, respectively. This choice results for 
example in 1200 imaginary time points and 2800 real 
time points for the lowest temperature run at /3 = 120 
(eV)“^ (T = 96 K), implying a total numerical grid of 
size 6800 x 6800 for a full double-time Keldysh contour 
Green function or self-energy. We checked that the nu¬ 
merical results were sufficiently converged as a function 
of time step size by adding runs with larger step sizes 
and extrapolating to zero step. Importantly, we used 
a 3rd order Adams-Bashforth scheme for the numerical 
integration of the differential equations together with a 
6th order Gregory integration for the numerical integrals. 
The self-consistency cycle for the self-energy typically re¬ 
quired a maximum of 5 iterations at each time step, us¬ 
ing a standard predictor-corrector scheme®®. Runs for 
the electron-phonon simulations typically required 40,000 
GPU hours. 


III. RESULTS 



FIG. 2. Light-enhanced superconductivity. Dynamics 
during and after r = 100 fs (dark colors) and r = 3 fs (light 
colors) ramps for different initial equilibrium temperatures 
(7c ,0 ~ 135 K). Solid (dashed) lines show the results of the 
el-ph (BCS) model and arrows indicate the Hnal thermal equi¬ 
librium values Af. Dark (light) grey shaded rectangles indi¬ 
cate the ramp durations. 

The fundamental question we address is whether and 
on which time scales superconducting order can be en¬ 
hanced by the change of the hopping amplitude. Fig. 2 
shows the dynamics of A(t) at different initial temper¬ 
atures for two different ramp times. For all situations 
a drastic enhancement of the order parameter is found 
for sizeable values of Aq. The initial increase is much 
faster for the short ramp duration of 3 fs compared to 
the slower ramp duration of 100 fs. After the ramp, 
the order parameter continues to increase for both ramp 
durations before it slowly approaches the thermal value 
Af. Damped “Higgs” amplitude mode oscillations are 
observed in some cases, as discussed in detail in Ref. 64. 

Importantly, the achieved enhancement of supercon¬ 
ductivity at short and intermediate times depends cru¬ 
cially on the initial order parameter and the distance 
from Af. To demonstrate the systematics. Fig. 3(a) 
shows the fraction of order parameter change during 
the T = 100 fs ramp, (Aj-amp — Ao)/(Ay — Aq), where 
Aramp = A{t = t) . The dashed line highlights the ap- 
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FIG. 3. Initial state and ramp duration dependence. 

(a) Order parameter change during the 100 fs ramp, Aramp — 
Ao, relative to A/ — Ao (symbols). Here Aramp denotes the 
order parameter at the end of the ramp, i.e. at the time t = 
T = lOOfs. The change scales almost linearly with the initial 
value Ao (dashed line), (b) Dependence of final steady state 
value on ramp duration r within BCS theory for different 
temperatures. Arrows show the data points for 100 fs ramps 
corresponding to Fig. 2. Dashed line indicates tAq = h. 


generate an effective attractive interaction between the 
electrons. 

The strong dependence of the order parameter dynam¬ 
ics on the ramp duration is already visible within BCS 
theory. Whereas the initial increase is strongly accel¬ 
erated for shorter ramp durations, the reachable steady 
state value®^ increases with longer ramp durations, as 
shown in Fig. 3(b), due to the reduced heating of the 
electrons and to the fact that slower ramps are less effi¬ 
cient at breaking Cooper pairs. The time scale to which 
the ramp duration has to be compared is Tf/Ag, around 
which the most important increase of the steady value 
takes place. The saturation value for long ramp dura¬ 
tions can be far from the thermal value at the temper¬ 
ature of the phononic bath. This failure to reach the 
thermal value is expected due to the presence of conser¬ 
vation laws within the integrable BCS model^®“^°’"^°. The 
exception is at T = 0 (top panel in Fig. 3(b)), where the 
order parameter follows closely the ground state value in 
the tAq oo limit even within the BCS approximation. 

In contrast, the full el-ph model exhibits a very dis¬ 
tinct behavior at intermediate and long times. The heat 
created in the electrons during the ramp is transferred 
to the phononic bath. As a consequence, the order pa¬ 
rameter at long times reaches the expected thermal value 
independent of the ramp duration. Surprisingly, as seen 
in Fig. 2, the phonon dissipative channel even allows to 
stabilize the strong initial increase directly after a rapid 
ramp. Thus fast ramps are more favorable than slow 
ramps with respect to light-induced superconductivity. 


proximate linear dependence of the achieved change on 
Aq at small Ag. In other words, the smaller the initial 
order parameter, the longer it takes to enhance super¬ 
conductivity. 

In order to gain a deeper understanding of the differ¬ 
ent regimes of the dynamics, we compare our results to 
simulations of a BCS model (see Appendix A1 A) with 
parameters chosen to match Ag and Ay in Fig. 2. The 
BCS model contains the electronic dynamics at a mean- 
field level including the phononic action only as an ef¬ 
fective pairing interaction between the electrons. Thus, 
the BCS model is not expected to be able to reproduce 
the full dynamics of the el-ph model. The purpose of the 
comparison between BCS and full el-ph dynamics is to 
illustrate the importance of dissipation of energy into the 
phononic bath. 

The BCS dynamics captures the main features of the 
initial increase of the order parameter in agreement with 
the el-ph model. However, small deviations occur in par¬ 
ticular for fast ramps, and the BCS model completely 
fails to account for the full thermalization on longer time 
scales. The reasonable agreement at initial times demon¬ 
strates that the dynamics of the order parameter at initial 
times is dominated by the change of the coherence factor 
relating the bare electrons to the Bogoliubov quasiparti¬ 
cles. In this regime, the main role of the phonons is to 




FIG. 4. Thermalization at long times via el-ph cou¬ 
pling. (a) Deviations of order parameter from respective final 
thermal values for various temperatures and full el-ph cou¬ 
pling parameter set (i) on a logarithmic scale. The approach 
to the final thermal value is well described by exponential 
decays (red lines), (b) Same as in (a) but for reduced el-ph 
coupling parameter set (ii). The slope in the exponential de¬ 
cays is indeed smaller by a factor of 0.8 as expected from the 
ratios of the electron-phonon couplings. Hence thermalization 
takes longer for smaller el-ph coupling. 

Fig. 4 shows the exponential decay of the order param¬ 
eter deviation from the final thermal values at long times 
on a logarithmic scale. Red lines indicate the slopes cor¬ 
responding to exponentially decaying behavior. As ex¬ 
pected the long-time relaxation, which is enabled by el- 
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FIG. 5. Time- and momentum-resolved dynamics. Snapshots of the momentum-resolved deviations from hnal thermal 
equilibrium values of the normal and anomalous occupations at times as indicated above, shown in 1/8 Brillouin zone each 
(separated by black lines) for a case of small Ao (a-c) and intermediate Ao (g-i). Effective Bogoliubov dispersions Ek{t) at the 
corresponding times along a. kx = ky momentum cut, compared to Hnal thermal dispersions Ekj (dashed curves) for the small 
(d-f) and large (j-1) Aq. 


ph coupling in particular to the acoustic phonon branch, 
crucially depends on the el-ph coupling strength. There¬ 
fore we compare parameter sets at different tempera¬ 
tures for full (a) and reduced (b) el-ph coupling. In¬ 
deed, we observe that the slopes are the same for dif¬ 
ferent curves within a panel, which have the same el-ph 
coupling, whereas the slopes are steeper for full compared 
to reduced el-ph coupling. In fact, the ratio of extracted 
slopes is 0.8 and matches approximately the ratio of the 
bare coupling values g'^, despite the fact that the self- 
energies are computed self-consistently and thus contain 
also higher orders in . The results clearly demonstrate 
the importance of el-ph coupling for the effective ther- 
malization of the superconducting state at long times. 

In order to gain additional insight into the interplay of 
collective order parameter dynamics and single-particle 
scattering during and after a 100 fs ramp, we show in 
Fig. 5 snapshots of the momentum-resolved dynamics of 
normal and anomalous densities, as well as the quasipar¬ 
ticle dispersions at selected times. We plot differences 
from the final thermal state, taken as the state with the 
final hopping value at the equilibrium temperature, to 
demonstrate the relaxation towards this state. 

Initially the distributions show strong deviations from 
the final thermal distribution in a broad momentum 
range. Right after the ramp, the normal distribution far 
from the Fermi surface quickly approaches the thermal 
distribution in both cases. Large deviations remain in a 
narrow region close to the Fermi surface. In contrast, the 
relaxation dynamics of the anomalous densities is much 
slower. In particular, for the case of smaller Aq almost no 
change can be detected, whereas for the intermediate Aq 
the momentum region of large deviations shrinks faster. 
Thus, the time scales for the relaxation of the normal 
densities are much faster than for the anomalous ones. 

This is further supported by snapshots of effec¬ 


tive Bogoliubov quasiparticle dispersions Ek{t) = 
y/e{k, ty + A(t)2, compared with the thermal disper¬ 
sions Ekj = ^efikY -b Aj (Fig. 5(d-f) and (j-1)), where 

e{k,t) = Ze{k,t). Whereas for small Aq (upper panels), 
the deviation from E^j is pronounced close to the Fermi 
surface, the dispersion for larger Aq (lower panels) is al¬ 
most thermalized at 392 fs. 

In order to highlight the nonthermal character of the 
instantaneous distributions, we plot in Fig. 6 the actual 
distributions compared with final thermal distributions 
as well as a reference set at much lower temperature of 
46 K. Importantly, the out-of-equilibrium normal distri¬ 
bution rife correspond to a “colder” fictitious instanta¬ 
neous temperature than both the final thermal and 46 
K reference distributions. The anomalous distribution 
/fc, on the other hand, is too “warm” in all cases. This 
discrepancy demonstrates that an effective quasi-thermal 
discription of the non-equilibrium data as often used, 
e.g. in the two temperature model, is entirely inade¬ 
quate here. This finding stresses the importance of a 
proper non-equilibrium modeling of the light-enhanced 
superconductivity. 

To emphasize the importance of the phonon spectrum, 
we compare in Fig. 7(a) the order parameter dynamics 
with and without the low-energy acoustic phonon branch. 
The corresponding Eliashberg functions are shown in the 
inset to Fig. 7(a). Clearly, the system effectively reaches 
the equilibration stage in the presence of the acoustic 
branch, whereas it is stuck at a nonthermal stationary 
value of the order parameter in the absence of the acous¬ 
tic branch. As discussed before at long time scales in the 
presence of the acoustic phonons shows an exponential 
relaxation towards the thermal state with a time scale 
approximately proportional to g^. 

To reveal the underlying reason for this nonthermal 
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behavior in the abscence of acoustic phonons, we show 
in Fig. 7(b) the momentum-resolved normal and anoma¬ 
lous density deviations from the hnal thermal values for 
the optical phonon. Clearly, there is a narrow win¬ 
dow around the Fermi surface without allowed scatter¬ 
ing phase space. This window is set by the optical 
phonon frequency and the Bogoliubov quasiparticle dis¬ 
persion £'fc(392 fs). A particle-hole pair with energy in 
the range [—Oopt/2, flopt/2] cannot relax because the re¬ 
quired energy transfer is smaller than flopt ■ This provides 
a very intuitive explanation for the importance of acous¬ 
tic phonons. Note, that the neglect of electron-electron 
scattering in our model is not the reason why the system 
shows this nonthermal behavior. In principle, it is cor¬ 
rect that the thermalization of electrons amongst each 
other after excitation is facilitated by electron-electron 
interactions. However, this possible thermalization, in 
the absence of phonons, would be bound to occur at a 
higher effective temperature than the initial equilibrium 
temperature, simply by energy conservation in the closed 
electronic system. Thus, light-enhanced superconductiv¬ 
ity would not profit from thermalization via electron- 
electron scattering. 


FIG. 6. Cuts of momentum distributions. Momentum 
distributions n*, and fk along a diagonal cut = ky for the 
same data as in Fig. 5 for T = 133 K (panels (a) and (b)) and 
T = 130 K ((c) and (d)), at times as indicated above. The 
black curves are instantaneous, non-equilibrium data, the red 
lines are bnal thermal distributions. Blue lines show thermal 
reference data for a lower temperature of 46 K. 




FIG. 7. Importance of acoustic phonons for thermal¬ 
ization. (a) Time-evolution of the order parameter in the 
presence of a narrow optic phonon mode (red curve) and ad¬ 
ditionally a broad acoustic phonon branch (blue curve). Inset: 
Eliashberg functions for the respective cases shifted for clar¬ 
ity. (b) Momentum-resolved deviations from final thermal 
values of normal and anomalous occupations at 392 fs in the 
presence of the narrow optic phonon mode (without acoustic 
phonons), shown in 1/8 of the Brillouin zone each separated 
by the black line. The relaxation processes are suppressed in 
a momentum window around the Fermi surface, which is de¬ 
termined by an energy window Ek{392 fs) < Hopt/2 (dashed 
curves). 


IV. CONCLUSIONS AND OUTLOOK 

In conclusion, we have demonstrated that nonequilib¬ 
rium superconductivity can be enhanced on time scales 
reachable in pump-probe experiments with THz pump 
pulses. An enhanced electronic density of states around 
the Fermi level leads to strengthening of the effective 
pairing interaction, which dynamically enhances the su¬ 
perconducting order parameter during and after a ramp 
of the electronic hopping amplitude. The main features 
of the short-time dynamics are well described by a BCS 
model which means that the presence of the phonons 
mainly enters via the effective attractive interaction. In 
contrast, the presence of a phononic bath with a broad 
spectrum to which electrons can release energy is cru¬ 
cial to ensure that the thermal state is reached in which 
the superconducting order is fully enhanced to its ex¬ 
pected equilibrium value after the ramp. Intriguingly, 
the phononic bath also enables the stabilization of the en¬ 
hanced order parameter already for very fast ramps and 
thus opens an interesting route towards light-enhanced 
superconductivity on very short time scales. 

The strong dependence of the dynamical enhancement 
of superconductivity on the initial order parameter raises 
the question of how to induce superconductivity when 
starting above Tc- A proper description of order pa¬ 
rameter fluctuations, which trigger the symmetry break¬ 
ing when starting in the normal state, is crucial in or¬ 
der to address this question (see Ref. 27 and references 
therein). Similarly, the nonequilibrium self-consistent 
update®®’®®’®® of the pairing phonons is an interesting 
topic for future research. 
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Al. APPENDIX 

A. Time-dependent BCS equations 

As a simplified alternative to the full Migdal- 
Eliashberg theory, we also use time-dependent BCS- 
theory in order to describe the initial evolution of the 
system. The BCS Hamiltonian is given by 

n =^e(fc,t)4^Cfc^ 

k(7 

- 1^1 4'A'iC-kiCkf, (Al) 

h,h' 

where U is the effective attractive interaction between 
the electrons, which is mediated by the electron-phonon 
interaction in the el-ph model. 

A mean-field decoupling in the Cooper channel leads 
to 


k(7 

+ + h.c. + constants. (A2) 

k 

We define the normal and anomalous densities 

= (cLcfca). (A3) 

fk = (c_fe 4 _Cfc.|.) = fk+ ifk y (A4) 

with real ( ) and imaginary part ( ). The time-evolution 
equations for these densities are given 

dtfkit) = Mk,t)fk{t), (A5) 

dtfkit) = -Mf^yt)fkit) - ^( 0(1 - rikit) - n-kit)), 

(A6) 

dt^il - nk{t) - n_k{t)) = 2A(t)4'(t). (A7) 

The self-consistency condition is 

Ait) = -\U\Y,fUtl (AS) 

k 


where we have used that we chose the initial equilibrium 
solution to be real and given by 


In e(fc,0)tanh(^) 

2 (1 - -’ 


.4(0) = 


-Aotanh(^M) 


2£;fc(0) 


4(o) = o, 


t.nnhf ^ 


2£;fc(0) 


(A9) 

(AlO) 

(All) 

(A12) 


Ek{0) = ^A2+e(fc,0)2. (A13) 


Here Ek{0) is the well-known Bogoliubov quasiparticle 
dispersion that is obtained from the diagonalization of 
the mean-field BCS Hamiltonian, where Aq is the self- 
consistently determined initial equilibrium order param¬ 
eter. 

The full solutions to the BCS equations are used in 
Fig. 2(a) and (c) of the main text. An analytical short- 
time solution for these equations can be obtained in 
the limit where we keep the normal density constant, 
nk{t) ~ nfe(O), and ignore self-consistent feedback by 
keeping A(t) Ri Aq fixed. The short-time limit of the 
obtained solution is 


4(t)-4(0) 


_ -Aq tanh( ^^^^°^ ) e(fc, 0)4^0 - 4) ^3 
3Afc(0) Jqt 

+ O((e(fc,0)t)4). (A14) 


This shows that for very short times t < tw, smaller 
than the inverse of the initial electronic half bandwidth 
{tw = 0.658 fs for lF/2 = 4Jo = 1 eV), the change 
in the momentum-resolved order parameter within BCS 
theory scales cubically in time for all momenta. For 
times t > tw the short-time approximation breaks down 
and the self-consistent feedback from other momenta be¬ 
comes crucial for the dynamics. Nevertheless, the ex¬ 
tracted scaling with Aq shown in Fig. 2(b) of the main 
text shows that the initial order parameter sets the im¬ 
portant dynamical time scale that governs the early-time 
enhancement of superconductivity during the ramp. 


B. Intermediate time behavior 

At intermediate times we observe a quasi-linear be¬ 
havior of the time-evolution of the order parameter, if its 
value is still far from final thermal value. This regime is 
reached, if the ramp duration is short compared to Aq, 
i.e. if tAq < h. 

Using linear fits to the intermediate time A{t) (see 
Fig. Al(a)) 

Afit(t) = A(161 fs) -f a{t - 161 fs), (A15) 

we obtain effective rates of change of the order param¬ 
eter. In Fig. Al(b) we show the rate normalized by 
the remaining deviations from the final thermal value, 
i.e. the effective rate of change irhaKAf — A(161 fs)). 
It is an interesting observation that the rates scale al¬ 
most linearly with the instantaneous order parameter for 
small A(I6I fs). By contrast, the rates deviate from this 
linear behavior at larger A(161 fs). In the regime of 
deviation, the rates also depend on the el-ph coupling 
strength. This behavior is in agreement with the ob¬ 
served momentum-dependent relaxation discussed in the 
context of Fig. 3 in the main text, and again illustrates 
the interplay of slow order parameter evolution far from 
thermalization and fast single-particle scattering becom¬ 
ing relevant for the order parameter dynamics close to 
thermalization. 
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4(161 fs) [eV] 


FIG. Al. Intermediate time behavior, (a) Time evo¬ 
lution of the order parameter for ramp duration r = 100 
fs (grey shaded area). An approximately linear behavior is 
found for intermediate times (red shaded area). Here the 
temperature is varied as indicated, data are for the “1.0 g^’’ 
el-ph coupling parameters, (b) Rates of change a, obtained 
from fits to A(t) at intermediate times, scaled by remaining 
deviation A/ — A(161 fs), versus instantaneous order param¬ 
eter. Colored arrows point to the parameter sets at T = 133 
K (magenta) and 128 K (orange). The black line indicates 
the equality between the scaled rate and instantaneous order 
parameter. 


Note that the linear change of A(t) in the intermediate 
regime is an empirical observation in our full simulations 
for the electron-phonon system, and is for instance not 
matched by BCS results. It only holds for a certain tem¬ 
poral regime, and only for “fast ramps” with small initial 
Aq, where the system is relatively far from its thermal 
value after the ramp. Thus, we identify heuristically a 
time scale for the increase of the order parameter at inter¬ 
mediate times, which will be important for experimental 
realizations. 




